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A driven Ising model with friction due to magnetic correlations has recently been proposed by 
Kadau et al [Phys. Rev. Lett. 101, 137205 (2008)]. The non-equilibrium phase transition present 
in this system is investigated in detail using analytical methods as well as Monte Carlo simulations. 
In the limit of high driving velocities v the model shows mean field behavior due to dimensional 
reduction and can be solved exactly for various geometries. The simulations are performed with three 
different single spin flip rates: the common Metropolis and Glauber rates as well as a multiplicative 
rate. Due to the non-equilibrium nature of the model all rates lead to different critical temperatures 
at > 0, while the exact solution matches the multiplicative rate. Finally, the cross-over from Ising 
to mean field behavior as function of velocity and system size is analysed in one and two dimensions. 

PACS numbers: 05.50.+q, 68.35.Rh, 04.20.Jb 



1. INTRODUCTION 

Magnetic contributions to friction due to spin correla- 
tions have attracted increasing interest in recent years. 
One interesting aspect is the energy dissipation due to 
spin waves in magnetic force microscopy, where magnetic 
structures are investigated by moving a magnetic tip over 
a surface [Hill [3]. On the other hand, magnetic friction is 
also present in bulk magnetic systems which are in close 
proximity. In this context, Kadau et al. [4 proposed a 
simple model for magnetic friction mediated solely by 
spin degrees of freedom. In this model an Ising spin sys- 
tem is moved over a second spin system with constant 
velocity v along a boundary. This permanent pertur- 
bation drives the system to a steady state far away from 
equilibrium, leading to a permanent energy flow from the 
boundary to the heat bath. 

This problem can be analyzed for several different ge- 
ometries in one, two and three dimensions, as shown in 
Fig. [T] Besides the original problem of two half-infinite 
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Figure 1: (Color online) Overview of the geometries consid- 
ered in this work. The grey regions are the magnetic systems, 
while the green (dark) regions are the moving boundaries. 
The arrows indicate the motion of the subsystems. 



two dimensional systems moving along the one dimen- 
sional boundary, denoted 2db in the following, we will 
consider the homogeneous cases Id and 2d where all spins 
are at the boundary, as well as the experimentally rele- 
vant three dimensional case 3db. Additionally, we will 
extend the analysis to sheared systems in two [5l [H [7] 
and three [8j dimensions, denoted 1+ld and 2+ld. These 
systems are experimentally accessible within the frame- 
work of shear flow in binary liquid mixtures (for a review, 
see [9 ), though with conserved order parameter, while we 
deal with a non-conserved order parameter. 

This model has some similarities to the driven lattice 
gas (DLG) proposed by Katz et al. [10] (see [TT] for a 
review) , where a system is driven out of equilibrium by an 
applied field which favors the motion of particles in one 
direction. We will discuss these similarities throughout 
this work. 

The paper is organized as follows: In the first part we 
will introduce the model and geometries and present, in 
the second part, an exact solution of the model in the 
limit of high driving velocities v ^ oo, which will be 
checked numerically in the last part using Monte Carlo 
simulations. There we will also investigate the case of 
finite velocities v. 



11. MODEL 

Let us start with the simplest case denoted Id in Fig.[l] 
and consider two Ising chains with spin variables cr = ±1, 
nearest neighbor coupling K = pj {jS = l/k^T and we 
set /cb = 1) and L\\ sites each, interacting with boundary 
coupling i^b = /3Jh and moving along each other with 
relative velocity v. In the Monte Carlo simulation the 
upper system is moved v times by one lattice constant 
ao with respect to the lower system during each random 
sequential Monte Carlo sweep (MCS). As one MCS cor- 
responds to a typical spin relaxation time to = O(10~^s) 
[12] and ao = O(10~^^m), the velocity v is given in nat- 
ural units ao/to = (9(lcm/s) (we will set ao = to = 1 in 
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Figure 2: (Color online) Sketch of geometry Id after A = 2 
moves. Spin ctq,/ interacts with spin a 1^1-^2 with coupling Jb 
(green (gray) lines), while all other couplings are J (black 
lines) . 



the following). 

To simplify the implementation, instead of moving the 
upper part of the lattice with respect to the lower part we 
reorder the couplings at the boundary with time. This 
procedure is analogous to the Lees-Edwards or mov- 
ing boundary condition in molecular dynamics simula- 
tions of fluids [13] and leads to a system as shown in 
Fig. [2] Assuming periodic boundary conditions (PBC) 
ak^i = crfe,^ mod L|| in the parallel direction, the time de- 
pendent Hamiltonian reads 



1 L|| 

o'k,iO'k,i-\-i- 

k=0 1=1 



"^b ^cro,^cri,Z+A(t) (1) 
^=1 



with the time dependent displacement 
A{t) vt. 



(2) 



The second geometry considered in this work is the 
2db case shown in Fig. |3) which already was investigated 
by Kadau et al. [4 . Here we have a square lattice with 
L\\ X sites and periodic boundary conditions in both 
directions, i.e., cFk.i = cr/e mod l^,/ mod l,, • Note that es- 
pecially o-L^,i = o-Q^i. The Hamiltonian of this system 
becomes 



k=l 1=1 



iO'k,i+i ^ K^^k(^k,i(^k+i,i+/\k{t) 

(3) 
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Figure 3: (Color online) Sketch of geometry 2db after A = 2 
moves. 



with A/c(t) = and K^^^ = K for all rows except row 
/c = 0, where the couplings to row /c = 1 are shifted 
with constant velocity Ao(t) = A(t) = vt. The coupling 
^±,0 = across the boundary is allowed to be different 
from K. For v = ^ and Jb = J = 1 this system simplifies 
to the 2d Ising model in equilibrium, which was solved 
exactly by Onsager [14 and shows a continuous phase 
transition at 



log(l + \/2) 



2.2691853.. 



(4) 



Note that both systems are translationally invariant in || 
direction under the transformation / ^ / + 1 and obey 
reflection symmetry at the boundary under k ^ 1 — k. 



III. EXACT SOLUTION AT HIGH VELOCITIES 

In Ref. [4^ it was shown that for high velocities v ^ 1 
the properties of the 2db system become independent of 
V. This can be understood as follows: In the limit v ^ 00 
the interaction i^bC'"o,^cri,z+A(t) across the driven bound- 
ary becomes uncorrelated, as, in the Monte Carlo simu- 
lations, at large v the spin cri^/+A(t) is different in every 
trial step and can, for simplicity, be a randomly chosen 
spin cTi^rnd from row 1. Note that this simplification was 
checked within the simulations and indeed gave the same 
results, enabling us to perform simulations at v = 00. 
Thus the boundary coupling can be replaced by the ac- 
tion of a fluctuating boundary field /i, e.g.. 



(5) 



with stochastic variables fiki = ±1 {k = 0,1) under the 
constraint {fiki) = {(^ki) = ^b, where mb denotes the 
magnetization at the driven boundary. Here we used the 
translation symmetry (aki) = rrik and the reflection sym- 
metry at the boundary, rrik = ^i-/c- In Fig. [4] this map- 
ping of the driven system onto a system with fluctuating 
boundary fields is illustrated for the Id case. The next 
step will be to map the fluctuating fields onto static fields 
by integrating out the degrees of freedom fij^i. 



A. Ising model in a fluctuating field 

Consider a general Ising model with arbitrary cou- 
plings Kij in a static external field h^^^ and additional 
fluctuating fields of strength ki (note the factor /3 in all 
field quantities) 



KijCJiCjj - ^(/^r^ + h^J^i)(Ji (6) 



where the = ±1 are stochastic variables at site i with 
given average 



{lii) = rrii. 



(7) 
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Figure 4: (Color online) Mapping of the Id driven system, 
shown for A = 2, on two disconnected Id systems with fluc- 
tuating fields. 



Comparing Eqs. (10) and (12), we conclude that under 
the condition 

tanh bi = tanh ki (13) 
the partition function Z can be expressed in terms of ^eq? 



= n 



cosh ki 
cosh bj 



(14) 



Eq. jlaj 



To summarize, the coupling with strength ki to fluc- 
tuating fields = ±1 with given average (/i^) = can 
be written as coupling to static effective fields bi with 
strength given by Eq. (13). In the next section we will 
use this mapping to exactly solve the driven Ising model 
for high velocities v ^ oo. 



As this condition is given a priori^ averages containing 
fii can be calculated using the trace formula 



(8) 



with the probability distribution Pi{/J.i) = (1 + /iimi)/2, 
as then 



as assumed. With the decomposition 



rrii 



(9) 



the degrees of freedom fi in the partition function Z can 
be traced out, 

Z = Tr^^e-^^^ = Tr^e-^^^Tr^ JJe^^^^^^ 

i 

= Yi cosh ki Tr^ e"^'^° R + ^^^^ 



B. Application to the driven Ising model 



The general condition Eq. (13) for the effective static 
fields bi simplifies for the systems considered in this work: 
As all boundary spins are equivalent, rrii = mb, with 
coupling ki = leading to a uniform effective field 
h\) = artanh(mb tanhi^Tb) at the boundary, as we assume 
no additional external fields, h^^^ = 0. Inserting this 
into the equilibrium expression for the boundary magne- 
tization mb,eq(K, ^b) = SlnZeq/9/ib, wc cud with the 
self-consistence condition 



mb,eq[K, artanh(mb tanhi^Tb)] = mb 



(15) 



for the non-equilibrium order parameter mb- 

As 1 = 9mb,eq/Smb|mb=o cit criticality, we obtain 
a very useful connection between the reduced zero 
field boundary susceptibility of the equilibrium model 
Xheqi^) = 9mb,eq/S^bUb=o the Critical tempera- 
ture Tc of the driven system by expanding Eq. (15) to 
first order around mb = 0, namely 



X^U^c) tanhi^b,c = l. 



(16) 



In the following we will apply these results to the one and 
two dimensional model introduced in Section HIl 



where we used the fact that = ±1. 

On the other hand, the Hamiltonian of the equilibrium 
Ising model without fluctuating fields in a static field hi 
can be written as 

P^eq = - ^ Kij(7i(7j - ^ hidi = pHo - ^ ^di (11) 

i<j i i 

with Ho from Eq. if we let bi = hi - hf\ The 
partition function of this model clearly fulfills 

Zed = e-^^^^ = Tr^ e'^^^ ]J e^^^^ 

= JJ cosh bi Tr^ e"^'^° 11 + ^^^^ • 



C. Id case 

The effective Hamiltonian of the system Id in a fluc- 
tuating field reads 



pn 



E 

1=1 



Kaicji+i + {h^^' + K^lii)cji. (17) 



Applying the self-consistence condition Eq. (15) to the 



well known expression for the equilibrium magnetization 
of the Id Ising model [cf. [15] 



rrieciiK, h) 



sinh h 



\Je-^^ + sinh^ h 



(18) 
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we obtain the zero field magnetization of the Id driven 
system in the ordered phase for velocity v ^ oo, 



cosh 2Ki, - coth 2K 



cosh 2i^b — 1 
with critical temperature fulfilling 



e^^^ tanhi^b,c = 1, 



(19) 



(20) 



as Xeq{K) = e^^ in this case. Interestingly, Eq. (19) 



is equal to the spontaneous surface magnetization of the 
2d equilibrium Ising model [T6^, Chapter VI, Eq. 5.20] 
if we identify K and with the couplings || and ± 
to the surface, and consequently has the identical criti- 
cal temperature T^. For the special case K = i^b this 
gives the well known value from Eq. Q. However, we 
regard this equality as coincidence without deeper mean- 
ing, as Eq. (19) is solution of a simple quadratic equation 
with small integer coefficients when written in the natural 
variables. Nevertheless, we checked this identity in the 
2d case and found that we do not get the surface magne- 
tization of the 3d system by the same procedure, as the 
critical temperature is Tc ~ 4.058 (Eq. J60|) instead of 
the correct value = 4.511424(53) ^mi- 

To calculate other quantities we use the transfer matrix 
(TM) formulation: the TM of the Id equilibrium Ising 
model reads [cf. [T5] 



^K+h g-K 



-K 



(21) 



and the partition function of a periodic system with L\\ 
spins can be expressed as 



^eq — Tr Teq . 



(22) 



Using Eq. (14) and the conditions Eq. (13) we can write 
Z = TrT^ii (23) 
with the TM (we set h^^^ — from now on) 



^ coshi^b ^ 
cosh h ^'^ 

which can be written as 



tanh h=m tanh K]^ 



rr. u / e^(l + sm?/^) e ^ costp 
T = coshAb _K I Kn • /^ 
* e ^ cosy; [1 — simp) 



usmg 



sin ip = m tanh i^b • 



(24) 



(25) 



(26) 



The angle decreases from i/; = 7r/2atT = 0toi/' = 
at T > Tc. The eigenvalues of T fulfill 



T\t^) = \^\t^) 



(27) 



and are given by 



Ao,i 



T <Tr 



coshi^b(e^±e"'') T>T, 



(28) 



where Aq denotes the larger eigenvalue dominant in the 
thermodynamic limit. Note that in this limit the analog 
to the free energy density 



logAo = -(J + Jb) 



(29) 



of the driven system is simply a constant in the ordered 
phase T < Tc [35 . Nevertheless, we can calculate physi- 
cal quantities within this TM notation using expectation 
values, as the whole information of the half-infinite sys- 
tem is contained in the normalized eigenvectors 



\to) = 



\ti) = 



- sm ^ 
cos ( 



(30) 



with cos2(/) = m. Using the normalized TM T = T/Aq 
and the Pauli matrix M = diag(l,— 1), the magnetiza- 
tion, Eq. (19), can be expressed as 



m={to\M\to), (31) 
while the correlation function in || direction becomes 

^11 H 



{cricri+n) - {o'i){cri^n) 
{to\MT^M\to) - {to\M\tof 
A?Ao-"(to|M|tl)^ 



as T- = A-|t^)(t^|. We get the result 



m )e 
tanh" K 



T<T, 

T>Tr 



leading to the inverse correlation length 

Ao f 2ifb T < Tr 



^n' = log 



Ai 



log coth JsT T>Tc 



(32) 



(33) 



(34) 



Note that does not diverge at the critical point, a 
feature which would lead to a correlation length exponent 
u = 0. In Section IV we will argue that in finite systems 



the spin fiuctuations are not only mediated by the spins 
(Ji but also by the self consistent field m which fiuctuates 
at finite Ly , an effect which vanishes in the exact solution, 
as oo. 

From the nearest neighbor correlation function we can 
calculate the internal energy ey = — J((j^cr^+i) in || direc- 
tion 



sinh 2K sinh K^, 
— JtanhiiT 



T<T, 
T>T, 



(35) 
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Figure 5: (Color online) I ntern al energy e||(T), Eq. (35), and 
C||(T), Eq. (36), of the Id driven system at 

2'\ 



and the 

dashed lines are results for the one dimensional Ising model 
in equilibrium. 



Figure 6: (Color online) Spin flip probability A, Eq. (43), and 
energy dissipation rate P, Eq. (47), versus reduced temper- 
ature T/Tc for the Id system at v ^ oo, together with MC 
data for L\\ =2^^. 



as well as the specific heat cm = de\\/dT in || direction 



for the multiplicative rate p^^i^E) 



2K^ 



sinh^ K 



(cothi^b-l) T<Tc 



T>T, 



(36) 



cosh K 

On the other hand, the internal energy in ± direction is 
simply given by 

e±_ = —Ji^m^ (37) 

as the related spins are uncorrelated. 

Now we turn to dynamical properties of this system 
under a concrete MC Glauber dynamics (see Sec. |IVB| 
for details) and calculate the spin flip acceptance rate 
^ = (Pflip) and the energy dissipation rate P = dE/dt: 
Let (OCCr) denote the probability of picking a spin a 
with direction C =T7i and left and right neighbors a£^r 
with direction Q^r- These probabilities can be calculated 
using the matrices P| = diag(l,0) and P| = diag(0, 1), 
e.g., (TTi) = (tolPTTPjTPilto). As the third coupling 
partner /i of spin a, with direction ^/i, is uncorrelated 
at infinite velocity, the probability of a particular spin 
configuration becomes 

(CKCr)(C^) = (to|Pc.tPctPcJto){to|PcJio). (38) 

The spin flip probability of a given configuration is 
Pflip(A£:), with AE = AEi + AE2 = 2J(j{(ti + a^) + 
2 Jbcr/i, and A becomes the sum over all 2^ possible cases 

^= E Pflip(A^)(OCCr)(CM), (39) 
which can be written as 

= E ^c(e"'''^(C) + (C)) (40) 

C=T,i 



IV B 



Eq. (69), 



Pfl.p(A£;i)p^.p(A£;2) introduced in Sec. 
using the abbreviation 



= E P?iip(A^l){CKCr) 

C«,Cr = T,i 

= e-'^'iCCC) + 2e-2^{CCC) + (CCC). (41) 



Note that the two terms in Eq. (40) are equal and the 
acceptance rate is independent of spin ^ because m is 
stationary. The resulting acceptance rate becomes 



A= I 



' cosh(i^ + i^b) - sinh(i^ - i^b) 
4g2(K+Kb) sinh K cosh^ K sinh i^b 

^ e-^b coshi^b(l - tanhi^)2 
which simplifies for J = Jb to 

r e-'^^coth2i^ 



T<Tc 
T>T, 



(42) 



sinh 2K 



^ coshK 



T<T, 



T>Te 



(43) 



The calculation of the energy dissipation rate P per 
spin is very similar to the acceptance rate A (Eq. (|40|) 
and gives 

P = -2Jb E Xc(e-^^^(C)-(C)). (44) 

C=T,i 

Furthermore, P j A can be calculated for arbitrary dimen- 
sions and geometries, as it is solely a property of the 
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fluctuating field. We find 



P 

A 



-2Jb 



X<.(e-2^.{C)-{C» 



2Jb(m 



Ec=T.i^c(e-'^^(C> + (C» 

^^-(C)-(C) 

2 ^ l)tanhirb 



tanh i^Tb 



(45) 



For the magnetization Eq. ( 19 ) of the Id system this gives 
2J.e-- ^^^^ 

(46) 



P 



tanh i^b 
2Jbtanhi^b T > Tc 



which, multiphed with A from Eq. (71) and for J = J\)^ 
becomes 



2e-^^ coth2K 
tanh sinh 2K 

2e-^^ tanhK 
coshi^ 



T<T, 



T>Tr 



(47) 



These results are shown in Fig. [6) together with data 
from MC simulations. Note that these results are only 
valid for the multiplicative rate p^-p from Eq. (69). 

Finally we list the critical exponents for the Ta driven 
system at ^ oo to be 



1 

2' 



7 = 1, a = 0. 



(48) 



The behavior of this system at finite velocities v will be 
discussed in Sec. IIVI 



D. 2db case 

The 2db case can be solved exactly using the expression 
for the equilibrium surface magnetization mb,eq(^,^b) 
of the 2d Ising model in a static surface field h\) ob- 
tained by McCoy and Wu [16, Chapter VI, Eq. 5.1], 
with z = tanhi^ and ?/b = tanh/ib- The integral rep- 
resentation given in their work can be further evaluated 
and written in closed form, the results are given in Ap- 
pendix |Aj Eq. (A2). If we again use Eq. ( p!5| ) and set 
Vh = ^b^b, with Zb = tanhi^b, we can calculate the 
non-equilibrium boundary magnetization m\){z,z\)) nu- 
merically as solution of the self-consistence condition 



mb,eq(^,mb2;b) mb, 



(49) 



which is shown for J = 1 and several values of Jb in 
Fig. [T] The critical temperature Tc of the system can be 
evaluated from the reduced zero field boundary suscepti- 
bility Xb,eq(^), Eq. (|A3|, to give 
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Figure 7: (Color online) Boundary magnetization mb(T), 
Eq. (49), of the 2db system for J = 1 and several values of 
Jb- For Jb = the mb reduces to the surface magnetization 
of the 2d equilibrium Ising model, Eq. (A5). 



for the case Jb = J = 1 using Xb°eq(^c)^b,c = 1 (Eq. ([l6|). 



As the critical temperature Tc, Eq. (50), is larger than 
the equilibrium critical temperature Tc,eq = 2.26918 . . ., 
the driven boundary induces a surface phase transition 
where only the driven surface has long range order above 
Tc,eq. The velocity dependence of this transition and the 
resulting phase diagram is discussed in more detail in 
Section [IVI 



2.6614725655752 . , 



(50) 



E. l+ld sheared case 

If the motion of the lattice described by Eq. ([s]) is not 
restricted to one row but applied to the whole system 
we get a system with uniform shear. Then all A/g(t) = 
A{t) = vt are equal, and we assume K_\_^k = K± to get 



0'k,lO'k,W + ^±C^fe,^C^/c+l,Z+A(t)• 
/e=l ^=1 

(51) 

Note that this system is translationally invariant in both 
directions, a fact that drastically simplifies the analysis 
of the critical behavior. 

Now we will investigate this system in the limit v ^ oo. 
Then each spin a^i interacts, as depicted in Fig. [8) with 
its neighbors cri^±i^i±^^t^ via fluctuating fields, while the 
interaction to the parallel neighbors (Jk,i±i remains un- 
changed. Thus the system decomposes into L± identical 
Id Ising models which again can be solved exactly: The 
coupling to two fluctuating fields /i^^i and /ii,2 with equal 
strength ki on each site can be traced similar to Eq. (10) 
to give 
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Figure 8: (Color online) Mapping of the 1+ld sheared sys- 
tem, shown for A = 2, on L± disconnected Id systems with 
fluctuating fields 



the order parameter of the sheared 1+ld system 

m(K||,Kx) 

1 - 2e^^\\ + 26^^11 Ve^^^ii - 1 + tanh^ 
tanh K± 

with critical temperature fulfilling 

2e^^ii '^ tanhi^xc = 1, 



(55) 



(56) 



which gives = 1/log (^^\/3 + vTz) = 3.46591... for 

J|| = J± = 1. A generalization of Eq. (52) from two to / 
fluctuating fields per spin is straightforward and leads to 
the general criticality condition 



X(°)(ifc)/tanhifb,c = l. 



(57) 



Although this geometry can be solved exactly at v = oo 
we expect the phase transition to be strongly anisotropic 
(see, e. g., [19 ) with two different correlation length ex- 
ponents i^ii > u±. In fact we found such behavior, with 
strong evidence for the exponents =3/2 and u± = 1/2, 
details on this will be published elsewhere [20^ . 



F. Other geometries 



Z 



= Tr^e-^^°Tr^]^]Je^^^^^-^^ 

i j = l 

= Trcr e~^^° Y[ [cosh ki + rriiai sinh kif 



n 



777/ ■ 

1 + ^iyr sinh2/c^ 



,(52) 



with 



Ci = - (1 - m,^ + (1 + m,^) cosh 2ki) . (53) 



Equating Eq. (52) with Eq. (12) we conclude that static 



fields bi can replace the fluctuating fields /i^j, with aver- 
age m,-, if 



tanh bj 



2m j sinh 2k j 



1 — m| + (1 + m|) cosh 2ki ' 



(54) 



The sheared system is translationally invariant in both 
directions, leading to homogeneous values = ki = 
and bi = h. Inserting Eq. (54) into Eq. (19) we get 



For two more cases we can derive highly accurate esti- 
mates for the critical temperature Tc of the driven system 
when V ^ oo^ namely the 2d Ising double layer [20 with 
Hamiltonian 

1 L|| L|| 

PH{t) = - \_^^klm{o-k,l,m+l + Cr/c,^+l,m) + 

k=0 1=1 m=l 

(58) 

and the experimentally relevant 2+ Id sheared case 

L± L|| L|| 

PH{t) = - [^\\^klm{o-k,l,m+l-^0-k,l+l,m)-^ 



k=l 1=1 m=l 



■ Kj_aklm(^k+l,l+A{t),- 



(59) 



both on simple cubic lattices: With Eq. (57) we can ex- 



press Tc using the high temperature series expansion for 
the reduced zero field susceptibility Xeq^(^) of the 2d 
Ising model, which was calculated to higher than 2000*^ 
order recently using a highly efficient polynomial time 
algorithm [21 . Using this extremely accurate result we 
find, for J = Jb = 1, the critical temperatures 

Tc = 4.058782423137980000987775040680 . . . (60) 

for the two 2d layers, and 

Tc = 5.264750414514743550598017203424 . . . (61) 
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for the 2+ Id sheared system with / = 2 analogous to 
the 1+ld sheared system. Note that due to the high 
accuracy of the series these values can be calculated to 
approximately 500 and 700 digits, respectively. 

Just for reference we also give the critical temperatures 
for two more cases: The experimentally relevant 3db case 
shown in Fig. [l] as well as the quite theoretical 3d case 
of two three dimensional systems in direct contact along 
the fourth dimension. In the 3db case we find = 4.8(1) 
using the 8*^ order high temperature series from Ref. [22", 
Tab. IV], while in the 3d case we obtain Tc = 5.983835(1) 
using the 32*^ order series from |23] . 

All these higher dimensional geometries are expected 
to show strongly anisotropic behavior with two (2d, 3db 
and 3d case) or possibly even three different exponents 
(2+ Id case), the reader is referred to Ref. [20 . 



IV. MONTE CARLO SIMULATIONS 




Figure 9: (Color online) Interactions of surface spin ai in the 
Id case 



fulfill 

^b,abs(T) cx (-r)^ (64a) 

Xb,abs(r) (X |r|-^ (64b) 

Cb(r) (X |r|-" (64c) 

with reduced temperature r = T/Tc — 1 and critical ex- 
ponents /3, 7 and a. But, before we present the results, 
we have to take a closer look at the used spin flip rates. 



Method 



B. An integrable algorithm 



We now describe the algorithms used to investigate 
the driven system: For finite velocities v we shift the 
boundary couplings by increasing A{t) from Eq. ^ after 
every N/v random sequential single spin flip attempts, 
where N denotes the total number of spins. Using 10^ — 
10^ MCS per temperature, we measured the following 
boundary properties: The boundary magnetization per 
spin and the energy per bond parallel to and across the 
boundary of a given configuration 



Mb 



^b = 



2L|| 



^ 1 ^11 

k=0 1=1 
J ^ ^" 
k=0 1=1 
Jb ^" 

E ^0,^cri,^+A(t) (62c) 
II 1=1 



2Lu 



as well as the corresponding bulk quantities. From these 
time dependent quantities we calculate the averages of 
the magnetization, reduced susceptibility. Binder cumu- 
lant, internal energy and specific heat at the boundary. 



^b,abs (|Mb|) 
Xb,abs = 2L|| {{M^) 



{\M^.\?) 



eb 

Cb 



1 



mi? 

{{El) - {E^f) 



(63a) 
(63b) 

(63c) 

(63d) 
(63e) 



Note that we have absorbed the factor /3 ^ into x- Near 
criticality these quantities show power law behavior and 



While equilibrium properties are most efficiently in- 
vestigated in Monte Carlo simulations using cluster al- 
gorithms, non-equilibrium systems have to be treated 
with random sequential single spin flip dynamics like the 
non-conserved Glauber dynamics [24] or the conserved 
Kawasaki dynamics [25 . The driven system is perma- 
nently under an external perturbation which drives it out 
of equilibrium, while the internal degrees of freedom are 
coupled to a heat bath in thermal equilibrium. From this 
coupling the spin flip probability pflip(A£^) of a given en- 
ergy change AE fulfills the detailed balance condition 



pmp{-AE) 

just like in the equilibrium case (for details, see [20 ). 



(65) 



The most common rates fulfilling Eq. (65) are the 



Metropolis rate [26 and the Glauber rate [24j . 

p^p(A^) = min(l,e-^^^). 



l + e/5^^' 



(66a) 
(66b) 



Using these rates in simulations of, e.g., the Id driven 
system, Eq. ([T]), it turns out that for all 'i; > the critical 
temperature Tc{v) depends on the used rate (see also 
Fig. 15): We find, for ^ oo and Jb = J = 1, the values 
= 1.910(2) and = 2.031(2) for the Metropolis 
and Glauber rate, respectively, while the exact solution 
Eq. ( 20 ) of the model presented in Section III gives = 



2.269.... Note that a similar dependency was recently 
found in the DLG by Kwak et al [27 . 

How can these discrepancies be understood? And can 
we construct a rate that matches the analytical treat- 
ment, i.e., has the same Tc? This is indeed possible: 
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Figure 10: (Color onlin e) Spin flip probabilities of the 
Metropolis rate E g. (|66a[ ) (dashed black line, circles), the 
Glauber rate Eq. ( |66b| ) (dotted blue line, squares), and the 
multiplicative rate Eq. (69) (red line, diamonds) for the Id 
system at criticality ( J = Jh — !)• 



Consider a microscopic change, i. e. a spin flip, of spin 
at the boundary (see Fig. [9]), with energy difference 



(Jj + 2J\)(Ji(Jr^ 

{j)' 



(67) 



AE2 



where the sum runs over the z neighbors of in the same 
subsystem (z = 2 in the Id case), while ar is from the 
other side of the moving boundary. The idea of the exact 
solution presented in the last section was to treat spin 
CT f as a fluctuating variable jii at site i with appropriate 
statistics. By contrast, correlations of different strength 
are introduced between the two subsystems by the rates 
Eq. (66), because the influence of spin cjr depends on 
the actual state of the z spins Gj . This can be seen most 
easily in the case of the Metropolis rate (Jb = J): if, e.g., 
(Ti = —(Jj then A£^i = —2zJ and = 1 independent 
of (Jr (note that AE2 = ±2 J), while in the parallel case 
{ai = (Jj) AEi = 2zJ and strongly depends on 



(see Fig. [To]). 

Fortunately, these rate-induced correlations can be 
completely eliminated by requiring that the flipping prob- 
ability is multiplicative^ 

Pflip(A^i + AE2) = pflip(A^i)pflip(AE2). (68) 

Clearly this condition is not satisfied for the rates 

2zJ + 2J) = 1, while 
(again we assume = J). 



e.g., 



in Eq. (66), 
p^p(-2^J)ip(2J) 

Instead, for simulations of driven systems we propose 
the rate 



^(AE-AE„i„) 



(69) 



which is uniquely defined by the detailed balance condi- 
tion, Eq. (65), and the multiplicity condition, Eq. (68) 



Figure 11: (Color online) Magnetization mabs(^), Eq. (64a), 
of the Id system at = 00 for several system sizes L\\ from 
Monte Carlo simulations, together with the exact solution 



Eq. (19) 



of AE at given geometry; this assures that p^-p(A£^) is 
maximal but never larger than one. For our example 
Eq. (67) we find A£^i^min = —2zJ and A^2,min = -2Jb 
to fumll Eq. ( [68| ). This new rate reproduces the calcu- 
lated critical temperatures in all considered geometries, 
e.g. T* = 2.269(1) for the Id case at v ^ 00. 

The resulting spin flip rates for the Id case at critical- 
ity are shown in Fig. [lO] Clearly, the multiplicative algo- 
rithm Eq. ( 69 ) has a smaller overall acceptance rate than 
Eqs. (66) and is thus slightly less efficient: A finite-size 
scaling analysis of the acceptance rate A = {pmp) at criti- 
cality in the Id case yields = 0.476(2), A^ = 0.366(2) 
and A* = 0.242(2) for the three algorithms, rendering 
this method roughly two times slower than the Metropo- 
lis algorithm. In fact, A^ = 3^2 - 4 = 0.24264... can be 
calculated exactly from Eq. (|43|. 

Note that the Metropolis and Glauber rates can be 
considered as many particle rates, as pflip depends on the 
many particle state of all coupling partners, while the 
multiplicative rate corresponds to a product of two par- 
ticle contributions. We believe that the dynamics gener- 
ated by the multiplicative rate is generally simpler than 
the one generated by Metropolis or Glauber rates, mak- 
ing an exact solution more feasible. Whether this differ- 
entiation only holds for the non-conserved Glauber dy- 
namics or also for the conserved Kawasaki dynamics is 
subject of future work. 

In the next two sections we will investigate finite-size 
effects in the Id case as well as the cross-over behavior 
at finite velocities v in the Id as well as in the 2db case. 
We first turn to the Id case. 



Id case 



|36j. The constant AE^i^i is the minimum possible value 



The exact solution presented in Section [III] was derived 
in the thermodynamic limit Ly 00, as we assumed a 
constant and non-fluctuating order parameter m in the 
self-consistence condition Eq. (15). This led to the result 
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Figure 12: (Color online) Finite-size scaling plot of the re- 
duced susceptibility Xabs(T), Eq. (64b), of the Id system for 
V = oo and system sizes Ly = 2^^, . . . , 2^^, together with the 
exact mean field finite-size scaling function (black line) from 
Ref. [28]. The correction factor C2 = 2.7. 



10.0 
5.0 

2.0 



1.0 



0.5 



0.2 



s 







■ 


v = A- 


♦ 


v=16 


▲ 


v = 64 


T 


v = 256 




v= 1024 



10"^ 10° 10^ lO'* 



10-2 10"^ 10° 10^ 10^ 10^ lO'* 



that the correlation length ^y, Eq. (34), remains finite 
at criticality. However, in a finite system the assump- 
tion m = const is not fulfilled and finite-size effects oc- 
cur, leading to a non trivial dependency of the physical 
quantities on system size. The fiuctuating order param- 
eter gives rise to additional correlations between spins at 
large distances not included in the exact solution. As 
the driven system shows mean field behavior, we can use 
the standard finite-size scaling theory for mean field sys- 
tems: Near criticality the correlation length parallel to 
the boundary fulfills ^||(t) (X |r|~^ii with critical expo- 
nent z^ii = 2/(ib, where (ib denotes the boundary dimen- 
sion. We have di^ = 1 in both the Id and the 2db case, 
leading to z/y = 2 in these cases. 

To illustrate these finite-size effects in the Id case, in 



Fig. 11 we show the magnetization mabs(^), Eq. (64a), 
as function of temperature for v = oo and several system 
sizes The exact solution, Eq. (19), is only approached 
in the limit L\\ 



oo. 



The finite-size scaling behavior is demonstrated exem- 
plarily for the susceptibility Xahs{T)^ Eq. (64b), which is 



shown in a finite-size scaling plot in Fig. 12 After resca- 



lation of the MC data in the usual way we indeed find 
the expected mean field exponents 7 = 1 and u\\ =2, 
furthermore the data falls onto the universal finite-size 
scaling function calculated in Ref. [28]. The same anal- 
ysis was performed for the magnetization mabs(^) and 
specific heat c(T), Eq. (64c), verifying the other two ex- 



ponents P = 1/2 and a = 0. 

In summary, the Id and the 2db systems with bound- 
ary dimension d]j = 1 have the standard mean field ex- 
ponents and fulfill the exponent relations 



2/3 + 7 = di^uii 



(70) 



We now turn to finite velocities v: Then the Id sys- 
tem always shows a cross-over from mean field to Ising 
behavior with increasing system size Lu. Only in the 



Figure 13: (Color online) Velocity dependent cross-over be- 
havior in the Id case. Shown is the rescaled width of the 
critical region 6rv^^'^ as function of the cross-over scaling 
variable L\\/v for several velocities v and several system sizes 
L|| = 2^, . . . , 2^^ (see text). The inset shows the correspond- 
ing cross-over of the effective correlation length exponent 
from = 1/2 (MF, dotted line) to = (Ising non- 
critical, dashed line). 



limit V ^ 00 the system undergoes a phase transition at 
finite temperatures. To investigate this velocity depen- 
dent cross-over, we measured the width dr of t he cr itical 
region by analysing the Binder cumulant Eq. (63c). Us- 



ing least square fits of the simulation data to the simple 
approximation 



[1 + tanh(f/Jr)] f < 



^b(T) 



1 
3 

1 

3 1 + f/Sr 



1 



(71) 



f > 



with f = T/Tc — 1 and fit parameters Tc and Sr, for 
several velocities v and system sizes Lu we determined 



St and plotted them in Fig. [T3| We find that the cross- 
over scaling variable is L\\/v in this case, while the y- 
axis has to be rescaled as Srv^^'^ to get the correct limit 
Ly^^" Sr = const with z^y = 2 in the limit v ^ 00. At 
finite V the width Sr stops shrinking at ^ 9'^, where 
denotes the cross-over system size, and only goes to 
zero for cx), indicating a sharp phase transition in 
this limit. The inset shows the effective exponent z^eff 
obtained from the logarithmic derivative, 

_i dlog Sr 



^eff 



SlogLii 



(72) 



whose value changes from = 1/2 (MF) to = 
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I 



0.10 
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0.00 



0.68 expi-LJ^^^c) 
MC data 



23 2^ 2^ 2^ 2^ 2^ 

Figure 14: (Color online) Influence of the system size L± on 
the critical point in the 2db case at = oo and L\\ = 256. The 
effective critical temperature Tc{L±) shifts to higher values if 
L± < 10^±,c (see text). 



(Ising) with growing system size. In the next section we 
will see that this behavior changes substantially in the 
2db case. 
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Figure 15: (Color online) Phase diagram of the 2db case. Be- 
low Tc,eq the two-dimensional bulk is ordered, while surface 
order is observed even above Tc,eq up to the velocity depen- 
dent phase boundary T^(v). The position of this boundary 
depends on the algorithm, the blue line holds for the multi- 
plicative rate, Eq. (69), while the thin red dotted line holds 
for the Metropolis rate, Eq. (66a). At fixed temperatures be- 
tween Tc,eq and Tc(v) a velocity driven phase transition is 
possible. The points are results from MC simulations. 



In the 2db case the moving boundary is coupled to a 
two-dimensional Ising model, which undergoes a phase 
transition at Tc^eq^ Eq. (|4|, independent of the velocity 
V. In addition, the moving boundary shows a boundary 
phase transition at temperature T^(v\ which grows with 
increasing v and eventually approaches the value given 
in Eq. (50) for ^ oo. As T^(v) > Tc^eq foi" all > we 
expect a boundary phase transition with paramagnetic 
bulk. Then the correlation length perpendicular to 
the boundary is finite at criticality and has the Ising value 



e±,c(^)=U(^c(^)), 



with [T6l 



Cq'm 



4K -2 log coth K 
logcothir-2ir 



T < Tc 
T > T, 



c,eq 
c,eq 



(73) 



(74) 



For that reason, in the finite-size scaling analysis it is 
sufficient for given v to simulate systems with varying 
length Lii while holding the height L± fixed at a value 
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oo 



L_\_ ^ ^±,c('^)- To illustrate this behavior, in Fig 
we show the effect of different values of L± for v = 
and L|| = 256. Only below Lj_ ^ 32 the system feels the 
finite width L^, resulting in a shift of the effective critical 
temperature Tc{L±) to higher values. The strength of 
the shift is proportional to the correlation function in 
± direction, {(Jq,i(^l^,i) oc exp(— Lx/^±,c)- The curves 
collapse for L± > 32 showing that a ratio Lj_/^j_^c ^ 10 
is sufficient, as ^±,c(oo) = 3.66323 ... in this case. 

We performed MC simulations and determined the 
critical temperatures for different velocities v by perform- 
ing a finite-size scaling analysis of the boundary suscep- 
tibility Xabs,b(^) and the boundary cumulant Uh{T). For 



the multiplicative algorithm, Eq. ( |69[ ), we used 400.000 
MC steps per temperature, while for the Metropolis algo- 
rithm 50.000 MC steps per temperature were used. The 
results are given in Tab. |l] and are compiled into a phase 
diagram of the 2db case shown in Fig. [15] An impor- 
tant aspect of this phase diagram is the possibility of a 
velocity driven non-equilibrium phase transition at fixed 
temperature (double arrow): While the system is para- 
magnetic at = and up to Vc{T) (thick blue line), the 
boundary shows long range order above that velocity. It 
would be interesting to see this transition in experiments, 
which could be performed in the corresponding geome- 
try 3db (see Fig. [T]), e.g., using two close rotating mag- 
nets slightly above the Curie temperature. The magnets 
should be isolating to avoid Eddy currents [4j. 
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2.644(3) 
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256 


2.654(2) 


2.44(2) 


1024 


2.659(2) 


2.45(2) 


oo 


2.661(1) 


2.45(2) 



Table I: Velocity dependent critical temperatures Tc{v) for the 
2db case, estimated using the multiplicative rate, Eq. (69), 



with 400.000 MC sweeps per temperature as well as using 
the Metropolis rate, Eq. (66a), with 50.000 MC sweeps per 
temperature. 
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Figure 16: (Color online) Velocity dependent cross-over be- 
havior in the 2db case. Shown is the rescaled width of the 
critical region Sr ^eci[Tc{v)] as function of the cross-over scal- 
ing variable L\\/^eci[Tc{v)] for several velocities v and different 
system sizes L\\ = 2^, . . . , 2"*^° (see text). The inset shows the 
corresponding cross-over of the effective correlation length ex- 
ponent i/eff from z^eff = 1 (Ising, dashed line) to z/eff = 2 (MF, 
dotted line). 



In the 2db case the cross-over scaling variable can be 
determined from the Tc{v) dependency discussed above. 
The correlation length ^eq at the critical point of the 
driven system, Tc{v)^ plays a key role: The system is 
Ising-like as long as correlations span the whole system 
in both directions || and ±, i.e. as long as the system size 
L|| is of the order of the bulk correlation length ^eq at 
the critical point Tc{v) of the driven system, leading to 
the cross-over scaling variable L||/^eq[^c('^)]- Again, the 
rescaling of the y-ax.is can be determined by requiring 
that a data collapse is obtained in the limit v ^ 0, leading 
to the expression (5r ^eq[^c('^)]7 as ^eq cancels in this case 
and we get the required condition L\\Sr = const ^ as ^eq 
^-z^eq [yi this limit, and z^eq = 1- 

The resulting cross-over scaling plot is shown in 



Fig. 16 For all finite v > the critical behavior changes 
from Ising to mean field at the cross-over system size 
^ 6^eq[^c('^)]- Below this value Sr shrinks accord- 
ing to 5r (X (Ising, dashed line), while above this 

]^ /2 

value (5r (X L|| ^ holds (MF, dotted line). As the shift 
exponent at small velocities, defined by 



T,{v)-T,{0)^v^ 



(75) 



-1/2 



is close to 1/2 we have, for small (x v ^ ^ v 

The shift exponent ^ = 1/2 has also been found in a field 
theoretical calculation of the 2+ld system [29^ . 



V. SUMMARY 

In this work we investigated a recently proposed driven 
Ising model with friction due to magnetic correlations. 
The non-equilibrium phase transition present in this sys- 
tem was investigated in detail using analytical methods 
and Monte Carlo simulations. In the far from equilib- 
rium limit of high driving velocities v ^ oo the model was 
solved exactly by integrating out the non-equilibrium de- 
grees of freedom. The resulting exact self-consistence 
equation was analysed for various geometries, leading in 
many cases to precise values of the critical temperature 
Tc of the non-equilibrium phase transition. In the limit 

V ^ oo the system always shows mean field behavior due 
to dimensional reduction, independent of geometry. In 
the simplest one dimensional case denoted Id a complete 
analysis of both equilibrium as well as non-equilibrium 
quantities has been presented. These exact results are 
another example of mean field critical behavior in an ex- 
actly solvable driven system, just as in the case of the 
DLG in a certain limit [30 . 

The analytic results were reproduced using a mul- 
tiplicative Monte Carlo rate originally introduced in 
[30] , which eliminates correlations due to many parti- 
cle dynamics introduced by the common Metropolis and 
Glauber rates. We claim that this algorithm is gener- 
ally favorable to the Metropolis and Glauber rates if an 
analytical treatment is considered. 

The finite-size effects naturally emerging in the simu- 
lations were analyzed using finite-size scaling techniques, 
a perfect agreement with exactly known universal finite- 
size scaling functions [28] were found. 

We analysed the critical behavior at finite velocities 
and studied the cross-over behavior from low to high 
velocities: We found that the Id system only has a phase 
transition in the thermodynamic limit for v = oo, while 
systems with finite v always become non-critical at the 
cross-over system size ~ 9v. On the contrary, the 
two-dimensional case 2db already has an Ising type phase 
transition at v = 0, which changes to mean field behavior 
for all finite > in the thermodynamic limit, at a cross- 
over length ^11^ ^ 6^eq[^c('^)]- In this sense, the velocity 

V is a relevant perturbation, always driving the system 
to a non-equilibrium state. 

The Id system changes from mean field to non-critical 
Ising universality, while the 2db case changes from Ising 
to mean field type with growing system size Ly. This 
somewhat puzzling fact can be understood in terms of 
the critical width Sr of the transition as follows: As 
in general Sr ex ^^^^^ at criticality, in the two dimen- 
sional Ising case Sr ex while in the mean field case 

— 1/2 

with one dimensional boundary (5r cx Ly ^ . Thirdly, 
cx Ly = const in the Id case at finite v. In the cross- 
over the actual critical width Sr is always governed by 
the largest contribution, and so at sufficiently large sys- 
tem size L\\ the contribution with smallest dominates 
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and determines the critical behavior. As consequence in 
both cases the effective inverse correlation length expo- 
nent changes from a larger value at small Ly to a 
smaller value at large Ly (1/2 ^ in the Id case, 1^1/2 
in the 2db case). 

Comparing the results to the driven lattice gas (DLG) 
fiO\ , we note that the DLG also shows a continuous non- 
equilibrium phase transition from an ordered to a disor- 
dered state at a critical temperature which grows with 
growing driving field. However, in the DLG the particle 
number is conserved, while we deal with a non-conserved 
magnetization. 

Finally some remarks on strongly anisotropic criti- 
cal behavior: The sheared system denoted 1+ld shows 
strongly anisotropic behavior at criticality and v ^ oo, 
with strong evidence for the correlation length exponents 
z^ll = 3/2 and u± = 1/2, details on this will be pub- 
lished elsewhere [20 . Remarkably, this is a rare case of 
an exactly solvable non-equilibrium system with strongly 
anisotropic critical behavior. 
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Appendix A: SURFACE MAGNETIZATION OF 
THE 2d ISING MODEL 



The equilibrium surface magnetization mb,eq of the 2d 
Ising model in a static surface field /ib obtained by McCoy 
and Wu [161 Chapter VI, Eq. 5.1] as well as the reduced 
zero field boundary susceptibility 



(0) 
^b,eq ~ 



dm 



b,eq 



dhv, 



(Al) 



can be written in closed form not present in the litera- 



is sometimes denoted Xii, and a high 



ture yet [3l] . Xh someuimes uenoueu Xll? 
temperature series expansion was derived up to 10*^ or- 
der in Ref. [32 and up to 23*^ order in Ref. [33 . As the 
expressions for anisotropic couplings and K± become 
way too complicated, we only give the results for the 
isotropic Ising model with K\\ = K± = K here: Using 
the definitions z = tanhiC, y = tanh/ib we find 



J 



^b,eq(^,^) 



Xb°eq(^) 



27T Attw 1 



n 



(1 



(l + 2i^-8^^) 



^ V 1 

z \ 

2^ K{l6w^) E(16^^ 



Attw 



by^ 

z 

1 

~ 4 



-,16^^ 



yl/2 _ y-1/2 

2(z-i -z) 



(A2) 



(A3) 



with the abbreviations and the complete elliptic integrals [37 of the 1®*, 2"^ 

and 3^^ kind, K(m), E(m) and n(n, m) Note that the 
^ _ ^v-*- ~ ^ ) (A4a) variable w is also used in high temperature series analysis 

(X^z^Y of the bulk zero field susceptibility [21]. For /ib = the 

\ — 2z — surface magnetization Eq. ( A2 ) reduces to the well known 

expression 

1 + 2^-^2 
b = — ^ — (A4c) 



l^z' 



9. ^r.^ /cosh2i^ - coth2ir 

?^ (A4d) "^■'"''^''V co,h2A--l - 
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